# This code reproduces Tables 2 and A10.

library(xtable)
library(dplyr)
library(mvtnorm)
library(estimatr)

set.seed(416)

load("data/mk_extract.Rdata")
# Table 2 -----------------------------------------------------------------
# while somewhat cumbersome, construct estimates by country due to differences in coding and design

## Brazil
bra_mod <- lm_robust(mtg_attend ~ -1 + as.factor(satis_trust_baseline), 
                     data = bra, cluster = clu, subset = ZZ == 1)
bra_top <- nrow(summary(bra_mod)$coef)
bra_dif <- bra_mod$coef[bra_top]-bra_mod$coef[1]
bra_dif_se <- sqrt(vcov(bra_mod)[1,1] + vcov(bra_mod)[bra_top, bra_top] + 
                     2 * vcov(bra_mod)[1, bra_top])

bra_ests <- c(coef(bra_mod)[c(1, bra_top)], bra_dif)
bra_ses <- c(summary(bra_mod)$coef[c(1, bra_top),2], bra_dif_se)

## Colombia
col_mod <- lm_robust(mtg_attend ~ -1 + as.factor(satis_trust_baseline), 
          data = col, cluster = clu, subset = ZZ == 1)
col_top <- nrow(summary(col_mod)$coef)
col_dif <- col_mod$coef[col_top]-col_mod$coef[1]
col_dif_se <- sqrt(vcov(col_mod)[1,1] + vcov(col_mod)[col_top, col_top] + 
                     2 * vcov(col_mod)[1, col_top])

col_ests <- c(coef(col_mod)[c(1, col_top)], col_dif)
col_ses <- c(summary(col_mod)$coef[c(1, col_top),2], col_dif_se)


## Liberia
lib_mod <- lm_robust(mtg_attend ~ satis_trust_baseline, 
                     data = lib, cluster = clu, subset = ZZ == 1)

lib_ests <- c(coef(lib_mod)[1] + coef(lib_mod)[2] * range(lib$satis_trust_baseline[lib$ZZ == 1]),
             diff(range(lib$satis_trust_baseline[lib$ZZ == 1])) * coef(lib_mod)[2])

lib_coefs <- rmvnorm(1e5, coef(lib_mod), vcov(lib_mod))
lib_ses <- c(sd(lib_coefs %*% c(1, min(lib$satis_trust_baseline[lib$ZZ == 1]))),
            sd(lib_coefs %*% c(1, max(lib$satis_trust_baseline[lib$ZZ == 1]))),
            sd(lib_coefs %*% c(1, max(lib$satis_trust_baseline[lib$ZZ == 1])) -
                 lib_coefs %*% c(1, min(lib$satis_trust_baseline[lib$ZZ == 1]))))


## Pakistan
pak_mod <- lm_robust(mtg_attend2 ~ -1 + as.factor(satis_trust_baseline), 
                     data = pak, cluster = clu, subset = ZZ == 1)
table(pak$satis_trust_baseline)
pak_top <- nrow(summary(pak_mod)$coef)
pak_dif <- pak_mod$coef[pak_top]-pak_mod$coef[1]
pak_dif_se <- sqrt(vcov(pak_mod)[1,1] + vcov(pak_mod)[pak_top, pak_top] + 
                     2 * vcov(pak_mod)[1, pak_top])

pak_ests <- c(coef(pak_mod)[c(1, pak_top)], pak_dif)
pak_ses <- c(summary(pak_mod)$coef[c(1, pak_top),2], pak_dif_se)


## Philippines--no panel survey
phl_mod <- lm_robust(mtg_attend ~ -1 + as.factor(satis_trust), 
                     data = phl, cluster = clu, subset = ZZ == 1)
phl_top <- nrow(summary(phl_mod)$coef)
phl_dif <- phl_mod$coef[phl_top]-phl_mod$coef[1]
phl_dif_se <- sqrt(vcov(phl_mod)[1,1] + vcov(phl_mod)[phl_top, phl_top] + 
                     2 * vcov(phl_mod)[1, phl_top])

phl_ests <- c(coef(phl_mod)[c(1, phl_top)], phl_dif)
phl_ses <- c(summary(phl_mod)$coef[c(1, phl_top),2], phl_dif_se)



## Uganda
uga_mod <- lm_robust(mtg_attend ~ -1 + as.factor(satis_trust), 
                     data = uga, cluster = clu, subset = ZZ == 1)
uga_top <- nrow(summary(uga_mod)$coef)
uga_dif <- uga_mod$coef[uga_top]-uga_mod$coef[1]
uga_dif_se <- sqrt(vcov(uga_mod)[1,1] + vcov(uga_mod)[uga_top, uga_top] + 
                     2 * vcov(uga_mod)[1, uga_top])

uga_ests <- c(coef(uga_mod)[c(1, uga_top)], uga_dif)
uga_ses <- c(summary(uga_mod)$coef[c(1, uga_top),2], uga_dif_se)


table_2 <- as.data.frame(rbind(lib_ests, lib_ses, uga_ests, uga_ses, col_ests, col_ses, phl_ests, phl_ses, pak_ests, pak_ses, bra_ests, bra_ses))
table_2 <- cbind(c("Liberia", "", "Uganda", "", "Colombia", "", "Philippines", "", "Pakistan", "", "Brazil", ""), table_2)

colnames(table_2) <- c("Country",
                       "P(Attend | Lowest Trust)", 
                       "P(Attend | Highest Trust)",
                       "Difference")

print(xtable(table_2, digits = 3), include.rownames = F, file = "results/Table_2.tex")



# Table A10 ----------------------------------------------------------------

# counts of respondents at endline and baseline
endline <- sapply(list(bra, col, lib, pak, phl, uga), FUN = function(x){nrow(x)})
baseline <- sapply(list(bra, col, pak, uga), FUN = function(x){sum(with(x, in_baseline), na.rm = T)})
baseline <- c(baseline[1:2], NA, baseline[3], NA, baseline[4])
# autocorrelation in trust in police
lib_clu <- lib %>% group_by(clu) %>%
  summarize(satis_trust = mean(satis_trust, na.rm = T),
            satis_trust_baseline = mean(satis_trust_baseline, na.rm = T))
autocor <- c(with(bra, cor(satis_trust, satis_trust_baseline, use = "complete.obs")),
             with(col, cor(satis_trust, satis_trust_baseline, use = "complete.obs")),
             with(lib_clu, cor(satis_trust, satis_trust_baseline, use = "complete.obs")),
             with(pak, cor(satis_trust, satis_trust_baseline, use = "complete.obs")),
             NA,
             with(uga, cor(satis_trust, satis_trust_baseline, use = "complete.obs")))
# itt on trust in police
get_itt1 <- function(data){
  out <- summary(lm_robust(satis_trust ~ ZZ, data = data, cluster = clu))$coef[2,c(1, 5:6)]
  return(paste0(round(out[1], 2), " [", round(out[2], 2), ", ", round(out[3], 2),"]"))
}
get_itt2 <- function(data){
  out <- summary(lm_robust(I(satis_trust-satis_trust_baseline) ~ ZZ, data = data, cluster = clu))$coef[2,c(1, 5:6)]
  return(paste0(round(out[1], 2), " [", round(out[2], 2), ", ", round(out[3], 2),"]"))
}

itt1 <- sapply(list(bra, col, lib, pak, phl, uga), FUN = get_itt1)
itt2 <- sapply(list(bra, col, lib, pak, uga), FUN = get_itt2)
itt2 <- c(itt2[1:4], NA, itt2[5])

table_a4 <- data.frame(Country = c("Brazil", "Colombia", "Liberia", "Pakistan", "Philippines", "Uganda"),
           `Survey Design` = c("Panel", "Panel", "Repeated cross-section", "Panel", "Endline only", "Panel"),
           `N in endline` = endline,
           `N in baseline` = baseline,
           `Autocorrelation` = autocor, 
           `ITT on Trust in Police` = itt1,
           `ITT on Change in Trust` = itt2)
print(xtable(table_a4), include.rownames = F, file = "results/Table_A10.tex")